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A new equilibrium thermodynamics analysis tool was built based on the CEA method 
using the OpenMDAO framework. The new tool provides forward and adjoint analytic 
derivatives for use with gradient based optimization algorithms. The new tool was validated 
against the original CEA code to ensure an accurate analysis and the analytic derivatives 
were validated against finite-difference approximations. Performance comparisons between 
analytic and finite difference methods showed a significant speed advantage for the analytic 
methods. To further test the new analysis tool, a sample optimization was performed to 
find the optimal air-fuel equivalence ratio, ¢, maximizing combustion temperature for a 
range of different pressures. Collectively, the results demonstrate the viability of the new 
tool to serve as the thermodynamic backbone for future work on a full propulsion modeling 
tool. 


Nomenclature 


quantity for the i'” element 

quantity for the 7” chemical species 

alternate subscript for the k‘”, chemical species 

ratio of specific heats 

lagrange multipliers 

residual function 

chemical potential energy 

air-fuel equivalence Ratio 

modified Lagrange multipliers 

density 

stoichiometric constant for the i*” element in the j*” species 
moles of an element summed across all species in a gas 
moles of an element summed across all species in a gas at the initial composition 
specific heat at constant pressure 

Lagrangian of the Gibbs free energy minimization 

Gibbs free energy 

specific enthalpy 
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° 


Al enthalpy of a species as a function of T 

ho input specified enthalpy 

n concentration of chemical species (kg-mol/kg-mixture) 
Ne total number of elements 

N, total number of chemical species 

P Pressure 

Ps Pressure at standard conditions (1.01325 Bar) 

R universal gas constant 

Ss entropy of a species as a function of T 

So input specified entropy 

T Temperature 

U vector of state variables in the chemical equilibrium solve 


I. Introduction 


A new generation of aircraft concepts have been proposed, exploiting propulsion airframe integration for 
higher efficiency. These concepts, such as over wing nacelledJ2] and boundary layer ingestion! couple the 
thermodynamic performance of the propulsion system with the aerodynamic performance of the airframe. 
This, in turn, necessitates a coupling of the aerodynamic and propulsion analysis tools for the aircraft design 
and motivates the application of Multidisciplinary Design Analysis and Optimization (MDAO) to navigate 
the potentially large design space. 

In order to apply MDAO, it is necessary to have a set of analysis tools well suited to optimization. 
This generally means that tools needs to be numerically stable, be capable of giving physically meaningful 
results even when given somewhat poor designs, and be reasonably computationally efficient. In cases where 
gradient based optimization will be employed, it is also important that analysis tools have smooth and 
differentiable formulations and it is highly beneficial if they provides analytic derivatives as well. 

The aerodynamics community has developed a set of Computational Fluid Dynamics (CFD) tools, such 
as Fun3d su25) and SumADE with adjoint analytic derivatives tailored for gradient based optimization. 
But the state-of-the-art propulsion analysis tool, the Numerical Propulsion System Simulation (NPSs) is 
not as well developed for optimization applications. NPSS is a flexible tool for modeling the thermodynamic 
cycles of propulsions systems and predicting their overall performance throughout a range of flight conditions. 
While it has been used in a number of integrated propulsion airframe studies, the methods used have been 
limited to lower fidelity aerodynamic models, loose coupling, and mostly gradient free optimizers. Felder 
et. al. applied NPSS to boundary layer ingestion and distributed propulsion concepts, though this did not 
directly couple the propulsion and airframe models! Allison et. al. have done extensive work integrating 
NPSS into the conceptual design for military aircraft, but have focused more on design space exploration 
than optimization S19] Geiselhart et. al. used NPSS as part of a low-boom supersonic transport design 
optimization, but noted numerical stability as a primary motivation for using gradient free optimization 
methods 

One of the major challenges with integrating NPSS into the gradient based optimization methods used for 
aerodynamic shape optimization is that it relies on finite-difference approximations to compute derivatives. 
The finite-difference method has two primary draw-backs: its computational cost grows linearly with the size 
of the design space and it can have numerical accuracy issues. Hendricks et. al. highlighted acute inaccuracy 
of finite difference gradients when optimizing a four stage power turbine, using meanline turbomachinery 
models built in NPSS#4!! To address the challenges of applying NPSS in MDAO applications, a new cycle 
analysis tool, named PyCycle, is being developed with the goal of maintaining the thermodynamic cycle 
analysis capabilities of NPSS but also providing analytic derivatives for MDAO applications. 

PyCycle is being built with the OpenMDAO framework #21 leveraging its built in solvers, data passing, 
and automatic multidisciplinary derivatives capability. Like NPSS, Pycycle will be composed of a set of 
propulsion elements (inlet, compressor, combustor, turbine, etc.) which can be linked together in a variety 
of ways to model a full propulsion system. Further examination of the underlying physics of an element 
reveals a finer structure. For example, calculations in the compressor element can be broken down into five 
steps: 


1. determine the exit pressure based on the input pressure ratio 
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solve for ideal exit thermodynamic conditions using the exit pressure and assuming constant entropy 
3. compute real exit enthalpy given efficiency 
4. solve for real exit thermodynamic conditions given exit enthalpy and pressure 


5. compute the shaft power required given the change in enthalpy across compressor 


Steps 2 and 4 are thermodynamic solves and steps 1, 3 and 5 are engineering relationships. The ther- 
modynamics solves compute properties of the working fluid (enthalpy, entropy, temperature, pressure, etc.) 
given any two thermodynamics states. All propulsion elements can be broken down into some set of thermo- 
dynamic solves separated by engineering calculations, although the specific combinations vary from element 
to element. This makes the thermodynamic solves a fundamental underlying requirement for any propulsion 
analysis code. This, before work can begin on the propulsion elements for PyCycle the underlying ther- 
modynamic solves must be implemented and validated. This paper presents the development of PyCycle’s 
thermodynamics module along with validation results for the analysis and analytic derivatives. 

There are a number of different well respected methods for computing accurate thermodynamics prop- 
erties for air and air-fuel mixtures, many with pre-existing implementations available. NPSS provides a set 
of different thermodynamics libraries (CEA, JANAF, ALLFUEL, GasTbl), which the user can select from 
at runtime. ALLFUEL and GasTbl are very computationally efficient, but are based on tabular thermody- 
namics data and offer lower accuracy. The CEA and JANAF libraries both solve for chemical equilibrium, 
but the JANAF library only considers a fixed set of species. The CEAZ library is by far the most general, 
allowing an arbitrary number of user defined chemical species. CEA’s generality made it the clear choice 
amongst the libraries available in NPSS. The major downside to using CEA was that its based on an Fortran 
implementation that can’t easily be augmented with analytic derivatives. One other option was considered, 
based on Cantera!#4! an object oriented chemical kinetics library that also provides equilibrium calculations. 
Chin et. al. used Cantera in a PyCycle prototype to predict performance for a two stage compression system, 
but found it to be very slow#5] Cantera computed equilibrium by running a dynamic simulation until the 
time derivatives damp out, which can take a larger number of internal iterations. That time-stepping method 
would also make it significantly more expensive to compute analytic derivatives with. Since neither CEA or 
Cantera were usable directly, a new thermodynamics implementation was built based on CEA scheme for 
minimization of Gibbs free energy. 

The rest of the paper is organized as follows: a summary of the Gibbs free energy equilibrium equations 
are presented in Section [A] A short discussion of the specifics of the numerical solver used to converge the 
chemical equilibrium solve is presented in Sectiorf{C] The analysis validation methodology and results are given 
in Section|IT]| In Section [I[V|numerical validation of the derivatives is presented along with some performance 
characterization of analytic vs finite difference derivatives. Finally, results from optimizations to find the 
maximum combustion temperature are presented in Section |V| that demonstrate the speed improvements 
achieved via the use of analytic derivatives. Future work will build on this study to develop the PyCycle 
propulsion system elements themselves and model a full thermodynamic cycle. 


Il. Chemical Equilibrium Equations 


A. Thermodynamic Properties Prediction 


CEA uses a two-step process to find the thermodynamic state of a gas. First it solves a system of nonlinear 
equations to minimize Gibbs free energy which gives the equilibrium composition for the gas at the prescribed 
state. Second an additional set of equations are solved, using the equilibrium composition, to compute the 
thermodynamic state of the gas. These equations are documented in detail by Gordon and McBride in their 
seminal paper on CEA but the methods are fundamental to this work!!! and the key chemical equilibrium 
equations are reproduced here for completeness. We refer the reader to the original CEA publications for 
details on the second step, computing the thermodynamic properties from the converged solution. 


B. Gibbs Free Energy Minimization 


The full thermodynamic state of a real gas can be defined by any two independent physical state variables 
(i.e. T, P, p, S, h). The CEA method provides calculations for three specific combinations: (1) temperature 
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and pressure-TP, (2) enthalpy and pressure-hP, (3) entropy and pressure-SP. The TP solve is the most 
fundamental. The hP and SP use the same equations as the TP solve augmented with T as a new state 
variable and an additional residual to drive the solution to the prescribed value of h or S. 


1. Temperature-Pressure Solve 


The TP solve is derived as follows: 


Gibbs free energy, g, is defined as 
Ns 


=> ( (4375) (1) 


j 


where N, is the number of chemical species, and n; and 1; are the concentrations (kg-mol/kg-mixture) and 
the chemical potential of the j‘” species respectively. The chemical potential is a function of temperature 
(T’), pressure (P), and composition (n) given by 


wi HRD) SAT) P 
RT RT ar In (=) tIn(n;) — In (do) ; (2) 


Both H° and Ss? are functions of T and a number of constants, co...cg, which are given for each species as 
input to the code. R is the universal gas constant. 
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Equilibrium composition is defined by the values of the concentration variables, n, that minimize Gibbs free 
energy subject to conservation of mass constraints. Mass is tracked on an elemental basis, given by 


Ns 
Rmass i =D (aijn;) =0 (5) 


where N, is the number of species in the mixture; b? is the amount of each element from the initial compo- 
sition; aj; is the stoichiometric constant for element 7 of species 7. Equation (5), defining the mass balance 
residuals, yields one constraint per element present in the mixture, giving N. — number of individual elements 
— constraints. To solve the mass constrained minimization of Gibbs free energy, we form a Lagrangian 


Nz Ne Ns 
G= S> (wynj) + ‘2 vi S°( (azjn4) (6) 
j=l i=1 j=l 


where 2; is the Lagrange multiplier for the i‘” element. We take derivatives with respect to the n and \ 
variables to build a system of non-linear equations whose solution will minimize the Lagrangian, as follows: 


Ne Ns 
6G = % Gx Ndi; )) 543 Y= (aijnj) — b? | 5A; = 0. (7) 
j=l i=1 \j=l 


Since 6G is linear in dn; and 6A;, we can split Equation into two sets of equations. We can get N, 
equations, 


Ne 
Rogibbs = bj + > (Aiasz) = 0, (8) 
=i 


as residuals representing Gibbs free energy. For convenience, because of the factors of 1/RT in Equation 
and Equation (4), we define an alternative Lagrange multiplier as follows: 





yeeS 


RT (9) 
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This gives an alternate form of Equation (8), 
jig 
Roibbs = oe — me (mez) = 0. (10) 
— 
Equations and yield a system with N, + Ne unknowns and N, + N; residual equations, which can 
then be converged with a numerical solver. 


2. Enthalpy-Pressure Solve 


The hP solve retains the state variables (n and 7) from the TP solve and the associated residuals from 
Equation and (10). In addition, it adds a new state variable T and a new residual to drive the computed 
enthalpy to match the the specified value of enthalpy, ho: 


Ns 


Rn = ho — S- (njH3(T)). (11) 


J=1 


3. Entropy-Pressure Solve 


Like the hP solve, the SP uses the same setup as TP with one additional state variable and residual. In this 
case, the new state variable is the prescribed entropy, Sp. The new residual drives the computed entropy to 
match the perscribed one: 


Ns ° Ns 
Rs =S,- ny (» ee —In (=) — In(n;) + in (>: a) ‘ (12) 
j=l 7 k 


Note that the pressure term is non-dimensionalized by pressure at standard atmosphere (P,=1.01325 Bar). 
The reference condition is necessary because Entropy is only really defined as a delta from a reference 
condition. 





C. CEA Modified Newton’s Method for Chemical Equilibrium 
1. Newton Convergence Scheme 


Gordon and Mcbride applied Newton’s method to converge the chemical equilibrium system 13] Netwon’s 
method, applied to the Gibbs free energy minimization, consists of successive solves of the linear system, 


dR 
Wau = -RW) (13) 


where U = (n,7) for a TP solve and U = (n,z,T) for hP and SP solves. AU is iteratively computed and 
applied until the residuals given by Equations (5), (10), (11), and converge to withing a given tolerance. 

Note that Equation involves computing ;, which via Equation (2), requires taking the natural log 
of n. In addition Equation involves taking the natural log of T. Hence it is important that the iterations 
do not push the values of n or T negative during convergence. A standard Newton’s method implementation 
does not provide a means of ensuring positive values for any of the state variables, and in fact for this 
problem it commonly produces negative values for n during convergence. To solve this issue, Gordon and 
Mcbride break up AU into two parts, for (n,T) and a variables. They then treat the (n, 7) component 
of the updates as Aln(n)/n and Aln(T)/T. This modifies the Newton update equation for the n and T 
variables as follows: 








Ney = Ne exp (ee) (14) 
Tia Teen () . (15) 


By using the exponential form, negative values from the newton solve are converted into multiplicative 

updates that are always positive, so assuming a positive initial guess the iterations can never push the values 

negative. This method has a secondary benefit of dramatically improving convergence of the overall model. 
The 7 update variables are treated normally, with a standard Newton update equation 


Trot = Te + Ar. (16) 
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2. Computing dR/dU 


dR/dU is computed analytically, with non-zero elements for the TP solve given by 








=i ae 
ARegivs 5 _ J Tm as (17) 
nk pe nt : 

AR givbs 5 
a (18) 

dR mass a 
———_ = aj;. 19 
dn; a J ( ) 


When doing an hP or SP solve, additional non-zero derivatives of Equation with respect to T are given 
by 


dRyitts 3 OHS O88 





— —. 2 
dT OT OT y) 
For an hP solve, the residual in Equation (11) contributes the following non zero terms: 

dRp 

ohh pH? 21 

et = 0H; (21) 
N 

dR, x aH? 
el ;{| 2 +H? }. 22 
aT 2m ( ar + (22) 


Similarly for an SP solve, the residual in Equation contributes the following non zero terms: 
dRg  S3(T) P at 
a =—R In P, In(n;) + in » nm, | —1, (23) 


N; 
dRg _@s(_ OH} 
a = (m ae (24) 


J=1 





All of these non-zero terms can be assembled into a matrix, which is then inverted with a direct method 
since the size of the matrix is at most (N, + Ne +1) x (Ns + Ne +1), with the largest term (N,) being on 
the order of 100 species even for very complex mixtures. 


III. Analysis Validation 


A. Analysis Validation Approach 


The thermodynamics module will serve as the foundation for PyCycle, so it is of paramount importance 
that it accurately model thermodynamic properties of air and air-fuel mixtures across a wide range of 
temperatures and pressures. The new code was validated against CEA predictions to ensure accuracy. A 
regular grid of temperatures and pressures, given in Table [I] were run in both CEA and PyCycle and the 
results compared at each test point. A total of 3600 different conditions were examined with temperatures 
ranging from 200 to 4800 degrees Rankine and pressures from 1 to 1500 psi. This grid of test points was 
run for four comparison cases, with equivalence ratios (¢) of 0, 0.016, 0.33, and 0.5 to provide a wide range 
of combustion conditions. Equivalence ratio (the ratio of actual fuel-to-air-ratio to the stoichiometric value) 
is a convenient way to express the amount of potential combustion in a manner that is independent of the 
specific fuel being used. 

To compute the validation data, PyCycle was setup with with air at the temperature and pressure 
conditions prescribed in Table[I]and then combusted at the the given ¢ with Jet-A, defined as a hydrocarbon 
fuel, C,2H23, with a stoichiometric fuel-air ratio of 0.06446 corresponding to ¢ = 1. 

For Pycycle the combustion was modeled as an adiabatic process by computing the overall enthalpy of the 
gas-fuel combination and holding it constant while solving for a chemical equilibrium composition. Since the 
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Table 1: Regular grid of temperature and pressure conditions used to validate PyCycle against 
CEA. A total of 3600 test points were run for each of the four ¢ settings 0, 0.016, 0.33, and 0.5 





low | high | step 


Temperature (degR) || 200 | 4800 | 200 
Pressure (PSI) 1 | 1500 | 10 
































total enthalpy calculation was needed for the PyCycle combustion model anyway, the corresponding CEA 
runs were setup with h and P using the enthalpy output from the corresponding PyCycle case. For both 
PyCycle and CEA a reduced set of 14 chemical species were considered, in order to reduce the number of 
trace species present in the converged equilibrium solves and to improve computational speed. The associated 
CEA input file can be found in Appendix[A] 

Two different types of validation were performed. First, the predicted chemical equilibrium compositions 
were compared to ensure that the proper amount of each chemical species was present at each test point. 
Next, the actual thermodynamic properties (P,T,p,h,5S,Cp,Cv,7y) were compared. For the composition 
validation the average discrepancy between PyCycle and CEA was 5.2710~° moles and the maximum error 
was .00011 moles. The average discrepancy in thermodynamic properties prediction was 0.03% and the 
maximum error was .52%. These results demonstrate a strong agreement between PyCycle and CEA over a 
very wide temperature and pressure range. 

One unavoidable source of discrepancy between PyCycle and CEA comes from implementation details 
related to what is done when composition variables, n;, trend toward 0. Taking logarithms and anti- 
logarithms of such small numbers is numerically challenging. CEA establishes a lower limit (user selectable 
with le-5 as the default) below which a chemical species is removed from the equilibrium solve all together 
and assumed to be 0. 

While this approach has some nice aspects — shrinking the size of the non-linear system and saving the 
cost of a logarithm evaluations on tiny quantities — it introduces some significant discontinuities into the 
solution. Discontinuities are problematic for optimization applications because they are a non-differentiable 
behavior. The result of this effect can be seen clearly in Figure [1] where the number of active species in 
the converged solution is shown for various different values of @ over a range of specified temperatures. 
The number of active species varies the most for ¢ = 0 (from 4 to 11) because at lower temperatures the 
composition of air stays relatively close to atmospheric, but NO, starts to form and dissociation starts to 
occur at the higher temperatures. Note that the data in Figure[I]is pressure-averaged over the entires range 
of pressured from Table[]] There was a very slight negative correlation between pressure and the number of 
active species in the data, but compared to the temperature effect it was negligible. If a larger set of species 
were considered, the pressure effect could be more pronounced. 

The discontinuities from the CEA approach were undesirable for PyCycle, so a a different approach was 
employed. Instead of removing species from the solution, a numerical limiting was introduced that made it 
increasingly difficult for the solver to drive composition variable values below the given lower limit — le-10 
by default. This approach gives a smooth and differentiable response, better suited for the purposes of 
optimization, but it does introduce a small source of discrepancy between CEA and PyCycle. Since the 
discrepancy deals with chemical species present only in minuscule amounts, this source of error was deemed 
more than acceptable. 


B. Chemical Equilibrium Validation 


The first validation check performed was to confirm that PyCycle returned the same overall chemical com- 
position as CEA over the prescribed test case range from Table[I] The comparison was done on a per species 
basis and measured with absolute differences. Air, even when combusted with ¢=0.5, is composed of over 
70% diatomic nitrogen, which is fairly non-reactive. This means that any other chemical species (e.g. CO2 or 
HO), while important to the thermodynamic properties, are only make up between 2% and 5% each. These 
small values necessitate the use of absolute error to measure the discrepancy between CEA and PyCycle and 
this requires knowledge of the actual value for n; in order to be meaningful. Figure |2)/shows the mean nj; 
value across the entire validation data set as a vertical blue line, as well as the ||n; cma — 2; Pycycle||2 as an 
orange curve. 
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Figure 1: The total number of active chemical species in CEA equilibrium solutions over the test 
case set from Table for a different values of ¢. The lower threshold on n;, before a species was 
removed, was set to 1le-10. 


In Figure |2| most of the discrepancies between CEA and PyCycle are at least one order of magnitude 
smaller than the mean concentrations, indicating very strong agreement between the two codes. For NH3 
and NOs the discrepancies are on the same order of magnitude as the mean value, but given the extremely 
small levels this did not cause significant changes in the thermodynamic properties. Some discrepancy is 
expected due to the different methods for handling small concentrations, as discussed above. 

In addition to considering the discrepancy of each species individually, Figure[3]examines how the norm of 
the absolute error, ||Rcka — NPycyclel||2, varies with T and P for the g=0.5 case. There is a clear correlation 
between increasing temperature — and to a much weaker extent pressure — and increasing discrepancy 
between the codes. Although the discrepancy grows from 3e-5 to 2e-4 moving from the lower left to the 
upper right corners, it is still very small. This trend also corresponds well with Figure [I] where CEA starts 
introducing small concentrations of new species as the temperature rises. Since CEA and PyCycle use 
different methods to handle these small amounts of species, some additional discrepancy in areas where this 
method is most active are expected. The results from Figure [2] and Figure [3] demonstrate that both codes 
are in strong agreement in terms of overall composition across the entire validation test set. The overall 
conclusion is that CEA and PyCycle both compute the same composition, within solver tolerance. 

The predicted composition vectors will not directly influence Pycycle propulsion elements. The elements 
will work with the predicted thermodynamic properties of the gas mixture. However its important that 
CEA and PyCycle agree on predicted composition because the thermodynamic properties are computed 
as a function of n. So the validation of predicted composition serves as a preliminary validation for the 
thermodynamics themselves. The validation of the thermodynamic properties is presented in the next 
section, but this data reinforces those results by ensuring that the calculations are fundamentally based on 
the same chemical compositions. 


C. Thermodynamic Properties Validation 


Comparison of the predicted thermodynamic state variables between CEA and PyCycle was done using 
relative measurements, since all values have magnitudes larger than 1. Figure [4] shows mean discrepancy, 
measured across the full test set, for entropy (.'), temperature, (J), enthalpy (h), density (p), pressure, (P), 
specific heat ratio (y), and specific heat at constant pressure (C;,). The mean error for all cases, across all 
properties, was .03%, and the maximum error was .52%. Note that since both h and P are set directly in 
the CEA run cases, from the output of the corresponding PyCycle cases, these properties have the lowest 
errors in Figure The other errors are larger, but still less than .1%. This demonstrates an extremely 
strong agreement between CEA and PyCycle with regard to calculations of the thermodynamic properties. 
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Figure 2: Mean difference between the equilibrium gas composition predicted by CEA and by 
PyCycle sampled over the test case set from Table [i} The dashed, red line indicate the value of the 
mean composition for each species to give a comparison between the magnitude of the error and 
the amount of species present. 


This ensures that PyCycle has a firm thermodynamic foundation upon which cycle modeling elements can 
be built. 


IV. Validating the Multidisciplinary Derivatives 


We define multidisciplinary derivatives as the total derivative taken of the objective of constraints, with 
respect to the design variables of the model. In other contexts these could also be called total derivatives or 
coupled derivatives. Admittedly, in this work, there is only a single engineering discipline, thermodynamics. 
Since the model is built up of components, you can think of the model as being built up from sub-disciplines 
(e.g. chemical equilibrium and gas properties calculations). In that sense, any derivative taken across an 
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Figure 3: ||ncBA — npyCycle||2 for ¢ = .5. This represents overall discrepancy between CEA and 
Pycycle at any point. There is clear trend for increasing — though still small — discrepancy as 


higher temperatures, and a much weaker correlation with increasing pressure. 
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Figure 4: Mean relative error between the thermodynamics states predicted by CEA and by PyCycle 
sampled over the test case set from table 
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OpenMDAO model as a multidisciplinary derivative. We use “multidisciplinary” to highlight the application, 
in OpenMDAO, to computing derivatives across a coupled multidisciplinary model. In contrast partial 
derivatives are taken across a single component (i.e. derivatives of component outputs with respect to its 
own input variables). OpenMDAO expects these values to be provided to it by components in the models, 
or it can approximate their values with finite-differences. Either way, these values are needed in order 
to compute the multidisciplinary derivatives. Given the partials the multidisciplinary derivatives are then 
computed by OpenMDAO by solving a linear system which accounts for all the data connections and coupling 
between components. 
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Figure 5: Finite-difference error for derivatives of h, S, p, Cp, y with respect to T and P, taken at 
standard day conditions. The data is presented in a Jacobian style format, with a single row for each 
response and column for each independent variable. Errors below 10~!° were achieved at various 
step sizes. Notably, minimum error was found at different step sizes for different combinations of 
variables. 


Since the motivation of this work was ultimately to build a tool well suited to gradient based optimization, 
derivative accuracy is just as important as analysis accuracy. This section presents the results of derivatives 
comparisons between the multidisciplinary derivatives and finite-difference approximations. Figure [5] shows 
log-log plots of finite-difference derivative error vs step-size for derivatives of h, S, p, Cp, y with respect to 
T and P, taken at standard day conditions. Each row represents derivatives of one response variable, and 
each column represents derivatives with respect to one independent variable. Strong agreement, with errors 
below 10~1!°, were achieved all of the derivatives at different step sizes using the central difference scheme. 

In Figure | the central difference scheme converges faster than the forward or backward schemes, and 
overall offers better accuracy over of a range of step sizes. This is expected, since central differencing is 
24 order accurate compared to 1% order accurate for the other two schemes. This accuracy is not free, 
as central differencing requires two function evaluations for each derivative. This compute cost shows up 
clearly in Figure [6} which presents the compute times for various methods to assemble the full Jacobian of 
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design variables and responses represented in Figure [5] The central difference scheme has twice the compute 
cost of the right and left difference schemes. Analytic derivatives are faster then the FD schemes, with 
the forward analytic method being the fastest. Each finite difference step requires the convergence of the 
full nonlinear model, however OpenMDAO computes the analytic derivatives using only linear solves via 
Martins and Hwang’s unified derivatives equations#® The additional cost of the adjoint vs the forward 
analytic derivatives method has to do with the number of linear solves for required for each one. There were 
2 independent variables, and 5 responses. So for forward mode, 2 linear solves were needed. For adjoint 
mode 5 were needed. With larger design spaces the adjoint method would be faster because the design 
variables far out number the number objective and constraints. 


Method Time (sec) 
Analytic Forward @ 0.0046 
Analytic Adjoint ® 0.0085 
FD Forward Diff. ® 0.0110 
FD Central Diff. O 0.0213 
FD Backward Diff. & 0.0110 


Figure 6: Wall times required to compute a full Jacobian for derivatives of h, S, p, Cy, y with 
respect to T and P. All analytic methods were faster than the finite-difference methods with 
forward analytic being the fastest because there were 2 design variables and 5 responses. All of the 
compute times were generated on an Apple Macbook laptop with a 2.6 Ghz dual-core Intel Core i5 
processor and 16 Gb of memory. 


V. Optimizations 


A. Equivalence Ratio Optimization 


We ran a series of unconstrained optimizations at a fixed gas temperature of 512°R, varying ¢ to maximize 
combustion temperature across different pressures from 15 psi to 1500 psi. Figure |7| shows the results of 
the optimizations. Somewhat counter-intuitively, the maximum temperature occurs for @ between 1.13 and 
1.08. Under the assumption of perfect combustion, the maximum temperature would occur at stoichiometric, 
@ = 1, with every molecule of diatomic oxygen being converted to water. However equilibrium calculations 
take into account disassociation effects, which simultaneously lowers the maximum achievable temperature 
and causes that temperature to occur at a richer gd The effect of dissociation becomes less severe at higher 
pressures because they tend to favor the creation of slightly larger molecules and the reaction more closely 
approximates ideal combustion. In Figure |7| this manifests in both the increasing maximum temperature 
and the decreasing optimal @ value as pressure increases. These results were run with both finite difference 
and analytic derivatives, but the results were essentially identical both in terms of overall execution time and 
final objective value. For this simple, one design variable, case analytic derivatives don’t provide significant 
speed improvements or accuracy advantage. 

These results establish that the best conditions for combustion, even accounting for equilibrium chemistry 
effects, are higher pressures. Therefore, one would expect that adding pressure as a design variable would 
yield the same result without the need to perform the parameter study across pressure. 


B. Equivalence Ratio & Pressure Optimization 


To verify this a second series of optimizations were run with both ¢ and P as design variables, again seeking 
to maximize combustion temperature to test this behavior. For this set of runs, the numerical solver tolerance 
for the chemical equilbrium was varied from le-8 to le-12 to test the sensitivity of the different gradient 
methods to depth of convergence. Both derivative methods were able to recover the expected result of 
@ = 1.128 and P = 1500 (psi). However, unlike the single design variable run, the analytic method provided 
a clear performance advantage over the finite-difference method. The finite-difference approximation was at 
least twice as expensive as the adjoint analytic method, as shown in Figure [8] One of the main challenges 
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Optimal ¢ & Temperature vs Pressure 


1.13 


1.08 


4370 


TCR) 


4200 


14 1500 
P (psi) 


Figure 7: Optimization of Equivalence Ratio, ¢, to find maximum combustion temperature at a 
range of different pressures. Note that maximum temperature is reached at a slightly rich ¢ due to 
dissociation affects, which become less pronounced at higher pressures. 


with using finite difference derivatives is the need to have tight tolerances set on any numerical solvers, 
such as the one used to converge chemical equilibrium. The data in Figure |8] quantifies this effect by 
comparing performance across a chemical equilibrium different solver tolerances. The analytic derivatives 
give a nearly constant computation time. However the finite-difference take between 5.5 seconds and 110 
seconds depending on solver tolerance and finite-difference step size. The wide range in compute times are 
a result of inaccurate derivatives forcing the optimizer to take more iterations to find the correct answer. 
Significantly, the behavior for the le-6 step size is not monotonic with respect to solver tolerance either. The 
worst performance came when using a step size of le-6 with tolerances between le-11 and le-10, but the 
times dropped back down for step sizes between le-10 and 1le-9. This indicates a high degree of numerical 
instability in the finite difference derivatives. For both step-sizes, at a solver tolerance of le-8 the compute 
times started to rise. Beyond that point, numerical noise prevented the optimizer from converging using the 
finite difference derivatives. This result clearly highlights the value of the analytic approach. Even for just 
2 design variables, the analytic derivatives provide both a significantly faster and more stable optimization. 


VI. Conclusion 


The growing need to incorporate propulsion cycle modeling into multidisciplinary optimizations for air- 
craft design, which rely on gradient based optimization methods and adjoint analytic derivatives, has moti- 
vated the development of a new propulsion cycle modeling tool, PyCycle. In order to develop PyCycle it was 
first necessary to build the core thermodynamics modeling capability, because that is a central module that 
enables building the cycle components such as a compressor, combustor, or turbine. This paper presents 
the development and validation of this thermodynamics module, based on an equilibrium chemistry method 
originally developed for the CEA code. Since of the calculations for propulsion system elements revolve 
around the thermodynamic calculations it was essential that both the analysis and derivative calculations of 
this module be heavily validated. The accuracy of the analysis was validated against the original CEA code 
and overall results agreed very well across a range of temperature and pressure conditions and at different 
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FD step = le-6 
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Figure 8: Optimization performance as a function of solver tolerance for the chemical equilibrium 
solve. The green line shows performance of the analytic derivatives method. The blue and orange 
lines show performance of the finite-difference derivatives for step-sizes of le-5 and 1le-6 respectively. 


equivalence ratios. The analytic derivatives were first validated against finite difference approximations and 
then further validated via usage in sample optimizations. 

From the derivative validation data and sample optimization runs, performance comparisons were made 
between finite difference and analytic methods of computing multidisciplinary derivatives. The analytic 
derivatives were significantly faster than the finite difference methods, even for a two variable optimization 
maximizing combustion temperature by varying pressure and equivalence ratio. In addition, the analytic 
derivatives made the performance of the optimization much less sensitive to the tolerance of the numer- 
ical solver for the chemical equilibrium calculations. The improvements in speed and accuracy clearly 
demonstrate the value of the analytic derivatives approach and provide solid motivation for the continued 
development of propulsion elements, analogous to those in NPSS, for Pycycle. 
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Appendix A CEA Species Definition 


SIX-CHARACTER REFERENCE-DATE CODES 


Letters* Reference Numbers 
g Glenn Research Center Month/year calculated 
j Chase:NIST-JANAF Thermochemical Month/year of table 


Tables. JPCRD 1998 Monograph 9. 


tpis  Gurvich:Thermodynamic Properties 


of Individual Substances Year of volume 
x TRC Thermodynamic Tables, 

Texas A&M Month/year of table 
bar Barin:Thermochemical Data of 

Pure Substances Year of volume 


coda CODATA Key Values for 
Thermodynamics Year of volume 


srd Standard Reference Data Year of JPCRD journal 


*NOTE: Upper-case letters indicate coefficients have not been 
recalculated since NASA TM-4513, 1993 ("old" polynomial form). 


ORDER OF SPECIES 


1) Gaseous products/reactants 
2) Condensed products/reactants 


! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 
! 3) Gaseous/condensed reactants only (starts with "Air") 
! 

! 

! 


thermo 
200.00 1000.00 6000.00 20000. 3/ 6/01 
Ar Ref-Elm. Spec: NSRDS-NBS 35 1971. 
3 g 3/98 AR 1.00 0.00 0.00 0.00 0.00 0 39.94800 0.000 
200.000 1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 
0.000000000D+00 0.000000000D+00 2.500000000D+00 0.000000000D+00 0.000000000D+00 
0.000000000D+00 0.000000000D+00 -7.453750000D+02 4.379674910D+00 
1000.000 6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 
2.010538475D+01-5 .992661070D-02 2.500069401D+00-3.992141160D-08 1.205272140D-11 
-1.819015576D-15 1.078576636D-19 -7.449939610D+02 4.379180110D+00 
6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 
-9.951265080D+08 6.458887260D+05-1.675894697D+02 2.319933363D-02-1.721080911D-06 
6 .531938460D-11-9.740147729D-16 -5.078300340D+06 1.465298484D+03 
CH30H Hf:TRC(6/87) w5030. Chen,1977. 
2g 7/00 C 1.00H 4.000 1.00 0.00 0.00 0 32.0418600 -200940.000 
200.000 1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 11435.277 
-2.416642886D+05 4.032147190D+03-2.046415436D+01 6.903698070D-02-7 .598932690D-05 
4.598208360D-08-1 . 158706744D-11 -4.433261170D+04 1.400142190D+02 
1000.000 6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 11435.277 
3.411570760D+06-1.345500201D+04 2.261407623D+01-2.141029179D-03 3.730050540D-07 
-3.498846390D-11 1.366073444D-15 5.636081560D+04-1 .277814279D+02 
CH4 TPIS 1991 v2 pti p44 pt2 p36. 
2 g 8/99 C 1.00H 4.00 0.00 0.00 0.00 0 16.04246 -74600.000 
200.000 1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10016.202 
-1.766850998D+05 2.786181020D+03-1.202577850D+01 3.917619290D-02-3.619054430D-05 
2.026853043D-08-4 . 976705490D-12 -2.331314360D+04 8.904322750D+01 
1000.000 6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10016.202 
3.730042760D+06-1 .383501485D+04 2.049107091D+01-1.961974759D-03 4.727313040D-07 
-3.728814690D-11 1.623737207D-15 7.532066910D+04-1 .219124889D+02 
ICH4 TPIS 1991 v2 pti p44 pt2 p36. 
2 g 8/99 IC 1.00IH 4.00 0.00 0.00 0.00 0 16.04246 -74600.000 
200.000 1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10016.202 
-1.766850998D+05 2.786181020D+03-1.202577850D+01 3.917619290D-02-3.619054430D-05 
2.026853043D-08-4 . 976705490D-12 -2.331314360D+04 8.904322750D+01 
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1000.000 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10016.202 


3.730042760D+06-1 .383501485D+04 2.049107091D+01-1.961974759D-03 4.727313040D-07 


-3.728814690D-11 1.623737207D-15 


C2H4 


2 


g 1/00 C 
200.000 


7.532066910D+04-1 .219124889D+02 
TRC w2600 4/88.JPCRD 1975 v4 p251. Chem Phys 1985 v98 pl. 

2.00H 4.00 0.00 0.00 0.00 0 28 .05316 52500.000 

1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10518.689 


-1.163605836D+05 2.554851510D+03-1.609746428D+01 6.625779320D-02-7 .885081860D-05 


5.125224820D-08-1.370340031D-11 


1000.000 


-6.176191070D+03 1.093338343D+02 
6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10518.689 


3.408763670D+06-1 .374847903D+04 2.365898074D+01-2.423804419D-03 4.431395660D-07 


-4.352683390D-11 1.775410633D-15 


Ic2 


H4 


8.820429380D+04-1.371278108D+02 
TRC w2600 4/88.JPCRD 1975 v4 p251. Chem Phys 1985 v98 pl. 


2g 1/00 IC 2.00IH 4.00 0.00 0.00 0.00 0 28 .05316 52500.000 


200.000 


1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10518.689 


-1.163605836D+05 2.554851510D+03-1.609746428D+01 6.625779320D-02-7 .885081860D-05 
5.125224820D-08-1.370340031D-11 -6.176191070D+03 1.093338343D+02 


1000.000 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10518.689 


3.408763670D+06-1 .374847903D+04 2.365898074D+01-2.423804419D-03 4.431395660D-07 
-4.352683390D-11 1.775410633D-15 8 .820429380D+04-1.371278108D+02 
C1i0H8,naphthale Naphthalene. Chen,JPCFD v8 n2 1979 p527-535. 


2 


me 
Z 


5. 
-1. 


8. 
il 
co2 

3 


HO2 
2 


g 8/93 C 
200.000 


10.00H 8.00 0.00 0.00 0.00 0 128.17052 150580.000 


1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 20713.076 


602920990D+05 6.237519290D+03-5.226157400D+01 2.397710630D-01-2.912272160D-04 


1000.000 


-854965710D-07-4.816685340D-11 -1.114753250D+04 2.972172708D+02 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 20713.076 


905864570D+06-3. 163144910D+04 7.030252590D+01-6 .018395960D-03 1.141923860D-06 
161432010D-10 4.891928210D-15 1.962512700D+05-4.347785692D+02 
IC10H8,naphthal Naphthalene. Chen, JPCFD v8 n2 1979 p527-535. 

2 g 8/93 IC 10.00IH 8.00 0.00 0.00 0.00 0 128.17052 150580.000 


200.000 


1000.000 


tpis79 C 
200.000 


1000.000 


1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 20713.076 


-602920990D+05 6.237519290D+03-5.226157400D+01 2.397710630D-01-2.912272160D-04 
.854965710D-07-4.816685340D-11 -1.114753250D+04 2.972172708D+02 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 20713.076 


-905864570D+06-3. 163144910D+04 7.030252590D+01-6.018395960D-03 1.141923860D-06 
-161432010D-10 4.891928210D-15 1.962512700D+05-4.347785692D+02 


TPIS 1979 v2 pt1 p25 pt2 p29. 
1.000 1.00 0.00 0.00 0.000 28.01010 -110535.196 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8671.104 


-489045326D+04-2 .922285939D+02 5.724527170D+00-8.176235030D-03 1.456903469D-05 
.087746302D-08 3.027941827D-12 -1.303131878D+04-7 .859241350D+00 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8671.104 


-619197250D+05-1.944704863D+03 5.916714180D+00-5.664282830D-04 1.398814540D-07 
.787680361D-11 9.620935570D-16 -2.466261084D+03-1 .387413108D+01 


6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8671.104 
868662960D+08-7 .500377840D+05 2.495474979D+02-3 .956351100D-02 3.297772080D-06 


g 9/99 C 
200.000 


1000.000 


-318409933D-10 1.998937948D-15 5.701421130D+06-2.060704786D+03 


TPIS 1991 v2 pt1 p27 pt2 p24. 
1.000 2.00 0.00 0.00 0.000 44.00950 = -393510.000 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 9365 .469 


-943650540D+04-6 .264116010D+02 5.301725240D+00 2.503813816D-03-2.127308728D-07 
-689988780D-10 2.849677801D-13 -4,.528198460D+04-7 .048279440D+00 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 9365.469 


.176962419D+05-1.788791477D+03 8.291523190D+00-9.223156780D-05 4.863676880D-09 
-891053312D-12 6.330036590D-16 -3.908350590D+04-2.652669281D+01 


6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 9365 .469 


g 6/97 H 
200.000 


1000.000 


-544423287D+09 1.016847056D+06-2.561405230D+02 3.369401080D-02-2.181184337D-06 
-991420840D-11-8 .842351500D-16 -8.043214510D+06 2.254177493D+03 


DO(H2) : JMolSpc,v33 1970 p147. NSRDS-NBS 3 SEC 6 1972. 
1.00 0.00 0.00 0.00 0.00 0 1.00794 217998 .828 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 


-000000000D+00 0.000000000D+00 2.500000000D+00 0.000000000D+00 0.000000000D+00 
-000000000D+00 0.000000000D+00 2.547370801D+04-4. 466828530D-01 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 


-078774250D+01-1.819354417D-01 2.500211817D+00-1.226512864D-07 3.732876330D-11 
-687744560D-15 3.410210197D-19 2.547486398D+04-4.481917770D-01 


6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 


g 6/97 IH 


200.000 


1000.000 


.173757694D+08-1.312035403D+05 3.399174200D+01-3.813999680D-03 2.432854837D-07 
-694275540D-12 9.644105630D-17 1.067638086D+06-2.742301051D+02 


DO(H2) : JMolSpc,v33 1970 p147. NSRDS-NBS 3 SEC 6 1972. 
1.00 0.00 0.00 0.00 0.00 0 1.00794 217998 .828 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 


-000000000D+00 0.000000000D+00 2.500000000D+00 0.000000000D+00 0.000000000D+00 
-000000000D+00 0.000000000D+00 2.547370801D+04-4 .. 466828530D-01 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 


.078774250D+01-1.819354417D-01 2.500211817D+00-1.226512864D-07 3.732876330D-11 
-687744560D-15 3.410210197D-19 2.547486398D+04-4.481917770D-01 


6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 


g 5/99 H 
200.000 


.173757694D+08-1.312035403D+05 3.399174200D+01-3.813999680D-03 2.432854837D-07 
-694275540D-12 9.644105630D-17 1.067638086D+06-2.742301051D+02 


Hills JCP 1984 v81 p4458. Jacox JPCRD 1988 v17 p303. 
1.000 2.00 0.00 0.00 0.00 0 33.00674 12552.000 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10002.162 


-7.598882540D+04 1.329383918D+03-4.677388240D+00 2.508308202D-02-3.006551588D-05 
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1.895600056D-08-4.828567390D-12 -5.809366430D+03 5.193602140D+01 


1000.000 6000 


-0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10002.162 


-1.810669724D+06 4.963192030D+03-1.039498992D+00 4.560148530D-03-1 .061859447D-06 
1.144567878D-10-4.763064160D-15 -3.194418740D+04 4.066850920D+01 


3 tpis78 H 2.00 


200.000 1000. 
4.078322810D+04-8. 
-1.202860160D-08 3. 
1000.000 6000. 
5.608123380D+05-8. 
5.936628250D-11-3. 
6000.000 20000. 
4.966716130D+08-3. 


-1.371809730D-11 1 
IH2 
3 tpis78 IH 2.00 


200.000 1000. 
4.078322810D+04-8. 
-1.202860160D-08 3. 
1000.000 6000. 
5.608123380D+05-8. 
5.936628250D-11-3. 
6000.000 20000. 
4.966716130D+08-3. 


-1.371809730D-11 1 
H20 
2 g 8/89 H- 2.00 
200.000 1000 
-3.947960830D+04 5 
4.955043490D-09-1 
1000.000 6000 
1.034972096D+06-2 
9.426468930D-11-4 
H202 
2g 6/99 H- 2.00 
200.000 1000 
-9.279533580D+04 1 
2.509255235D-08-6 
1000.000 6000 
1.489428027D+06-5 
6 .947265590D-12-4 


3 g 5/97 N 1.00 


200.000 1000. 
0.000000000D+00 0. 
0.000000000D+00 0. 

1000.000 6000. 


8.876501380D+04-1 


4.012657880D-11-2. 
6000.000 20000. 
5.475181050D+08-3. 


-1.098367709D-11 1 
NH3 
2 tpis89 N 1.00 


200.000 1000. 


-7.681226150D+04 1 


1.317385706D-08-3. 

1000.000 6000. 
2.452389535D+06-8. 
2.530923570D-12-3. 


3 tpis89 N 1.00 


200.000 1000. 


-1.143916503D+04 1 


-7.685111050D-09 2. 
1000.000 6000. 


2.239018716D+05-1 


-1.416076856D-11 9. 

6000.000 20000. 
-9.575303540D+08 5. 
2.912584076D-11-3. 


Nno2 
2g 4/99 N 1.00 
200.000 1000 
-5.642038780D+04 9 
9.145497730D-09-1 
1000.000 6000 
7.213001570D+05-3 
-7.611335900D-11 3 

NO3 
2 j12/64 N 1.00 
200.000 1000 


Ref-Elm. TPIS 1978 vi pti p107 pt2 p31. 

0.00 0.00 0.00 0.00 O 2.01588 0.000 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8468.102 
009185450D+02 8.214701670D+00-1.269714360D-02 1.753604930D-05 
368093160D-12 2.682484380D+03-3.043788660D+01 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8468.102 
371491340D+02 2.975363040D+00 1.252249930D-03-3.740718420D-07 
606995730D-15 5.339815850D+03-2.202764050D+00 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8468.102 
147448120D+05 7.983887500D+01-8.414504190D-03 4.753060440D-07 
-605374600D-16 2.488354660D+06-6 .695524190D+02 
Ref-Elm. TPIS 1978 vi pti p107 pt2 p31. 

0.00 0.00 0.00 0.00 0 2.01588 0.000 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8468.102 
009185450D+02 8.214701670D+00-1.269714360D-02 1.753604930D-05 
368093160D-12 2.682484380D+03-3 .043788660D+01 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8468.102 
371491340D+02 2.975363040D+00 1.252249930D-03-3.740718420D-07 
606995730D-15 5.339815850D+03-2.202764050D+00 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8468.102 
147448120D+05 7.983887500D+01-8.414504190D-03 4.753060440D-07 
-605374600D-16 2.488354660D+06-6 .695524190D+02 
CODATA 1989. JRNBS 1987 v92 p35. TRC tuv-25 10/88. 

Oo 1.00 0.00 0.00 0.00 0 18.01528 -241826.000 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 9904.092 
-755731020D+02 9.317826530D-01 7.222712860D-03-7 .342557370D-06 
-336933246D-12 -3.303974310D+04 1.724205775D+01 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 9904.092 
-412698562D+03 4.646110780D+00 2.291998307D-03-6 .836830480D-07 
-822380530D-15 -1.384286509D+04-7 .978148510D+00 
Cons:TPIS 1978 vi pti pi21.Hf:TPIS 1989 vi pti p127. 

O 2.00 0.00 0.00 0.00 0 34.01468 -135880.000 
-0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 11158.835 
-564748385D+03-5 .976460140D+00 3.270744520D-02-3 .932193260D-05 
-465045290D-12 -2.494004728D+04 5.877174180D+01 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 11158.835 
.170821780D+03 1.128204970D+01-8 .042397790D-05-1 .818383769D-08 
.827831900D-16 1.418251038D+04-4.650855660D+01 
Hf :CODATA1989. Spec:NSRDS-NBS 3 sec5 1975. 

0.00 0.00 0.00 0.00 0 14.00674 472680.000 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 
000000000D+00 2.500000000D+00 0.000000000D+00 0.000000000D+00 
000000000D+00 5.610463780D+04 4.193909320D+00 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 
-071231500D+02 2.362188287D+00 2.916720081D-04-1.729515100D-07 
677227571D-15 5.697351330D+04 4.865235790D+00 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6197 .428 
107574980D+05 6.916782740D+01-6.847988130D-03 3.827572400D-07 
-277986024D-16 2.550585618D+06-5 .848769710D+02 
TPIS 1989 p219. JRNBS 1968 v72A p207 for low T. 

H 3.00 0.00 0.00 0.00 0 17.03056 -45940.000 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10043.121 
-270951578D+03-3 .893229130D+00 2.145988418D-02-2.183766703D-05 
332322060D-12 -1.264886413D+04 4.366014940D+01 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10043.121 
040894240D+03 1.271346201D+01-3.980186580D-04 3.552502750D-08 
322700530D-16 4.386191960D+04-6 .462330250D+01 
TPIS 1978,1989 vi pti p326 pt2 p203. 

Oo 1.00 0.00 0.00 0.00 0 30.00614 91271.310 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 9179.110 
-536467592D+02 3.431468730D+00-2.668592368D-03 8.481399120D-06 
386797655D-12 9.098214410D+03 6.728727490D+00 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 9179.110 
-289651623D+03 5.433936030D+00-3.656034900D-04 9.880966450D-08 
380184620D-16 1.750317656D+04-8 .501667090D+00 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 9179.110 
912434480D+05-1.384566826D+02 1.694339403D-02-1.007351096D-06 
295109350D-16 -4.677501240D+06 1.242081218D+03 
Hf0,Cons: TPIS 1989 vi pti p332 pt2 p207. 

QO 2.00 0.00 0.00 0.00 0 46 .00554 34193.019 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10208.175 
-633085720D+02-2.434510974D+00 1.927760886D-02-1.874559328D-05 
-777647635D-12 -1.547925037D+03 4.067851340D+01 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10208.175 
-832615200D+03 1.113963285D+01-2.238062246D-03 6.547723430D-07 
-328361050D-15 2.502497403D+04-4.305129910D+01 
JPCRD 1998 Mono.9 p1607. 

QO 3.00 0.00 0.00 0.00 0 62.00494 71128.000 
0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10958.914 
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3.405398410D+04 2.266670652D+02-3.793081630D+00 4.170732700D-02-5.709913270D-05 
3.834158110D-08-1.021969284D-11 7.088112200D+03 4.273091810D+01 


1000.000 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 10958.914 


-3.943872710D+05-8 . 244263530D+02 1.061325843D+01-2.448749816D-04 5.406060320D-08 
-6.195466750D-12 2.870000149D-16 8.982011730D+03-3 .444666500D+01 


3 tpis78 N 
200.000 


Ref-Elm. TPIS 1978 vi pti p280 pt2 p207. 
2.00 0.00 0.00 0.00 0.000 28.01348 0.000 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8670. 104 


2.210371497D+04-3 .818461820D+02 6.082738360D+00-8.530914410D-03 1.384646189D-05 
-9.625793620D-09 2.519705809D-12 7 .108460860D+02-1 .076003316D+01 


1000.000 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8670.104 


5.877124060D+05-2.239249073D+03 6.066949220D+00-6.139685500D-04 1.491806679D-07 
-1.923105485D-11 1.061954386D-15 1.283210415D+04-1.586639599D+01 

6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8670.104 
8.310139160D+08-6 .420733540D+05 2.020264635D+02-3.065092046D-02 2.486903333D-06 
-9.705954110D-11 1.437538881D-15 4.938707040D+06-1 .672099736D+03 


3 g 5/97 0 
200.000 


DO(02):CJP v32 1954 p110. Spec:NSRDS-NBS 3 sect 1976. 
1.00 0.00 0.00 0.00 0.00 0 15.99940 249175.003 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6725.403 


-7.953611300D+03 1.607177787D+02 1.966226438D+00 1.013670310D-03-1.110415423D-06 
6.517507500D-10-1 .584779251D-13 2.840362437D+04 8.404241820D+00 


1000.000 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6725 .403 


2.619020262D+05-7 .298722030D+02 3.317177270D+00-4. 281334360D-04 1.036104594D-07 
-9.438304330D-12 2.725038297D-16 3.392428060D+04-6 .679585350D-01 
6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6725 .403 
1.779004264D+08-1 .082328257D+05 2.810778365D+01-2.975232262D-03 1.854997534D-07 
-5.796231540D-12 7.191720164D-17 8.890942630D+05-2. 181728151D+02 


3 g 5/97 I0 
200.000 


DO(02):CJP v32 1954 p110. Spec:NSRDS-NBS 3 sect 1976. 
1.00 0.00 0.00 0.00 0.00 0 15.99940 249175.003 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6725 .403 


-7.953611300D+03 1.607177787D+02 1.966226438D+00 1.013670310D-03-1.110415423D-06 
6.517507500D-10-1 .584779251D-13 2.840362437D+04 8.404241820D+00 


1000.000 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6725 .403 


2.619020262D+05-7 .298722030D+02 3.317177270D+00-4.281334360D-04 1.036104594D-07 
-9.438304330D-12 2.725038297D-16 3.392428060D+04-6 .679585350D-01 
6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 6725.403 
1.779004264D+08-1 .082328257D+05 2.810778365D+01-2.975232262D-03 1.854997534D-07 
-5.796231540D-12 7.191720164D-17 8.890942630D+05-2. 181728151D+02 


3 tpis78 0 
200.000 


TPIS 1978 pt1 p110 pt2 p37. 
1.00H 1.00 0.00 0.00 0.000 17.00734 39344.106 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8813. 106 


-1.998858990D+03 9.300136160D+01 3.050854229D+00 1.529529288D-03-3.157890998D-06 
3.315446180D-09-1 . 138762683D-12 3.239683480D+03 4.674110790D+00 


1000.000 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8813.106 


1.017393379D+06-2.509957276D+03 5.116547860D+00 1.305299930D-04-8 . 284322260D-08 
2.006475941D-11-1.556993656D-15 2.044487130D+04-1.101282337D+01 
6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8813.106 
2.847234193D+08-1 .859532612D+05 5.008240900D+01-5.142374980D-03 2.875536589D-07 
-8.228817960D-12 9.567229020D-17 1.468642377D+06-4.023555580D+02 


3 tpis89 0 
200.000 


TPIS 1989 vi pt1 p94 pt2 p9. 
2.00 0.00 0.00 0.00 0.000 31.99880 0.000 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8680. 104 


-3.425563420D+04 4.847000970D+02 1.119010961D+00 4.293889240D-03-6 .836300520D-07 
-2.023372700D-09 1.039040018D-12 -3.391454870D+03 1.849699470D+01 


1000.000 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8680.104 


-1.037939022D+06 2.344830282D+03 1.819732036D+00 1.267847582D-03-2.188067988D-07 
2.053719572D-11-8 . 193467050D-16 -1.689010929D+04 1.738716506D+01 

6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8680.104 
4.975294300D+08-2.866106874D+05 6.690352250D+01-6 .169959020D-03 3.016396027D-07 
-7.421416600D-12 7.278175770D-17 2.293554027D+06-5 .530621610D+02 


102 
3 tpis89 IO 
200.000 


TPIS 1989 vi pt1 p94 pt2 p9. 
2.00 0.00 0.00 0.00 0.000 31.99880 0.000 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8680. 104 


-3.425563420D+04 4.847000970D+02 1.119010961D+00 4.293889240D-03-6 .836300520D-07 
-2.023372700D-09 1.039040018D-12 -3.391454870D+03 1.849699470D+01 


1000.000 


6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8680.104 


-1.037939022D+06 2.344830282D+03 1.819732036D+00 1.267847582D-03-2.188067988D-07 
2.053719572D-11-8 . 193467050D-16 -1.689010929D+04 1.738716506D+01 

6000.000 20000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8680.104 
4.975294300D+08-2.866106874D+05 6.690352250D+01-6.169959020D-03 3.016396027D-07 
-7.421416600D-12 7.278175770D-17 2.293554027D+06-5 .530621610D+02 


H20(s) 
11 8/89 H 
200.000 


ICE. SANFORD GORDON, NASA TP-1906, 1982. 
2.000 1.00 0.00 0.00 0.00 1 18.01528 -299108.000 
273.150 7 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 0.000 


-3.89736288D+05 2.41502204D+03 6.09419068D+01 -8.47031137D-01 4.47768117D-03 
-1.06521862D-05 9.77189851D-09 0.00000000D+00 -5.39150616D+04 -2.07596481D+02 


H20(L) 
11 8/89 H 
273.150 


COX ET.AL.,CODATA KEY VALUES FOR THERMODYNAMICS.1989. 
2.000 1.00 0.00 0.00 0.00 1 18.01528 -285830.000 
600.000 7 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 13278.000 


8.72123781D+07 -1.39087511D+06 9.15729532D+03 -3.17596351D+01 6.13885076D-02 
-6.26798865D-05 2.64536349D-08 0.00000000D+00 6.43316523D+06 -4.90985319D+04 


END PRODUCTS 
Air 


Mole%:N2 78.084,02 20.9476,Ar .9365,C02 .0319.NASA TP1906 1982 
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2 g 9/95 N 1. 


200.000 


56160 .4959AR.00936C .00032 .00000 0 28.9651784 -125.530 
1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8649. 264 


1.009950160D+04-1.968275610D+02 5.009155110D+00-5.761013730D-03 1.066859930D-05 


-7.940297970D- 


1000.000 


09 2.185231910D-12 -1.767967310D+02-3 .921500990D+00 
6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8649. 264 


2.415214430D+05-1.257874600D+03 5.144558670D+00-2.138541790D-04 7.065227840D-08 


-1.071483490D-11 6.577800150D-16 6 .462263190D+03-8 . 147408670D+00 
Inert Air Mole%:N2 78.084,102 20.9476,Ar.9365,C02 .0319.NASA TP1906 1982 
2 g 9/95 N 1.561710.41959AR.00937C .00032 .00000 0 28 .9651784 -125.530 
200.000 1000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8649. 264 
1.009950160D+04-1.968275610D+02 5.009155110D+00-5.761013730D-03 1.066859930D-05 
-7.940297970D-09 2.185231910D-12 -1.767967310D+02-3 .921500990D+00 
1000.000 6000.0007 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 8649. 264 
2.415214430D+05-1.257874600D+03 5.144558670D+00-2.138541790D-04 7.065227840D-08 
-1.071483490D-11 6.577800150D-16 6 .462263190D+03-8 . 147408670D+00 
H2(L) Hydrogen. JANAF Prop.Ser.D,3/66. 
O jp3/66 H 2.00 0.00 0.00 0.00 0.00 1 2.01588 -9012.000 
20.270 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000 
JP_7 NASA C. Snyder - 9/17/2001. Hcomb = 18875.BTU/# 
0 g 5/95 C 1.00H 2.0044 0.00 0.00 0.00 1 14.03102 -4745 .000 
298.150 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000 
JP-7 NASA C. Snyder - 9/17/2001. Hcomb = 18875.BTU/# 
0 g 5/95 C 1.00H 2.0044 0.00 0.00 0.00 1 14.03102 -4745 .000 
298.150 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000 
Inert_JP_7 NASA C. Snyder - 9/17/2001. Hcomb = 18875.BTU/# 
0 g 5/95 IC 1.00I1H2.0044 0.00 0.00 0.00 1 14.03102 -4745 .000 
298.150 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000 
Inert_JP-7 NASA C. Snyder - 9/17/2001. Hcomb = 18875.BTU/# 
0 g 5/95 IC 1.00I1H2.0044 0.00 0.00 0.00 1 14.03102 -4745 .000 
298.150 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000 
02(L) Oxygen. NASA RP-1311 PtII 1996. 
0 g 4/95 0 2.00 0.00 0.00 0.00 0.00 1 31.99880 -12979.000 
90.170 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000 
Inert_02(L) Oxygen. NASA RP-1311 PtII 1996. 
0 g 4/95 I0 2.00 0.00 0.00 0.00 0.00 1 31.99880 -12979.000 
90.170 0.0000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000 
JP_4 Or RP-1(L). NASA RP-1311, Part II, 1996. Hcomb = 18640.BTU/# 
015/95 C 1.00H 1.9423 zl 13.96872 -22723.000 
298.15 
JP-4 Or RP-1(L). NASA RP-1311, Part II, 1996. Hcomb = 18640.BTU/# 
015/95 C 1.00H 1.9423 i 13.96872 -22723.000 
298.15 
JP-5 Or ASTMA1(L). NASA RP-1311, Part II, 1996. Hcomb = 18600.BTU/# 
015/95 C 1.00H 1.9185 z 13.94473 -22183.000 
298.15 
Jet-A(L) NASA TM-101475,1988. Hcomb=18500 BTU/#:NASA CR-72951,1971. 
11 2/96 C 12.00H 23.00 0.00 0.00 0.00 1 167 .31462 -303403.000 
220.000 550.000 7 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 0.000 
-4.218467340D+05-5 .576234840D+03 1.522094335D+02-8.610096140D-01 3.071640926D-03 
-4.702766120D-06 2.743009309D-09 0.000000000D+00-3 . 238535350D+04-6 . 780954740D+02 
Jet-A(g) NASA TM-101475,1988. NASA CR-72951,1971. 
21 2/96 C 12.00H 23.00 0.00 0.00 0.00 0 167 .31462 -249657 .000 
273.150 1000.000 7 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 0.000 
-6.068699280D+05 8.328264220D+03-4.312323550D+01 2.572391032D-01-2.629316827D-04 
1.644989491D-07-4.645336690D-11 0.000000000D+00-7 .606965040D+04 2.794307229D+02 
1000.000 5000.000 7 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 0.0 0.000 
1.541871660D+07-7 .433869020D+04 1.468645380D+02-1.297042936D-02 2.159140196D-06 
-1.887183642D-10 6.604559540D-15 0.000000000D+00 3.996323340D+05-9 .266674660D+02 
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